Introduction à la modélisation statistique bayésienne

Thierry Phénix & Ladislas Nalborczyk
18 juin 2019

Problème n°1

Peut-on prédire la taille d'un individu par la taille de ses parents ?

library(tidyverse)

d1 <- read.csv("Data/parents.csv")
glimpse(d1)

Observations: 40
Variables: 4
$ gender <fct> M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, …
$ height <dbl> 62.5, 64.6, 69.1, 73.9, 67.1, 64.4, 71.1, 71.0, 67.4, 69.…
$ mother <int> 66, 58, 66, 68, 64, 62, 66, 63, 64, 65, 64, 64, 62, 69, 6…
$ father <int> 70, 69, 64, 71, 68, 66, 74, 73, 62, 69, 67, 68, 72, 66, 7…

Solution problème n°1

La taille de la mère a l'air plus prédictive de la taille d'un individu, et ce d'autant plus si cet individu est une femme…

d1 %>%
    gather(parent, parent.height, 3:4) %>%
    ggplot(aes(x = parent.height, y = height, colour = parent, fill = parent) ) +
    geom_point() +
    stat_smooth(method = "lm", fullrange = TRUE) +
    facet_wrap(~gender) +
    theme_bw(base_size = 25)

plot of chunk unnamed-chunk-2

Solution problème n°1

On peut fitter plusieurs modèles avec brms, et les comparer en utilisant le WAIC.

library(brms)

m1.1 <- brm(
    height ~ 1 + gender,
    prior = c(
        prior(normal(70, 10), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 10), class = sigma)
    ),
    data = d1,
    iter = 4000, warmup = 1000, chains = 4, cores = 4)

m1.2 <- update(
    m1.1,
    newdata = d1,
    formula = height ~ 1 + gender + mother + father)

m1.3 <- update(
    m1.1,
    newdata = d1,
    formula = height ~ 1 + gender + mother + father + gender:mother)

m1.4 <- update(
    m1.1,
    newdata = d1,
    formula = height ~ 1 + gender + mother + father + gender:father)

Remarques :

On utilise la fonction update car on ne veut changer que lq formule
On introduit newdata car on ajoute les variables 'father' et 'mother'

Solution problème n°1

m1.3 %>%
    plot(
        pars = "^b_",
        combo = c("hist", "trace"), widths = c(1, 2.5),
        theme = theme_bw(base_size = 20)
        )

plot of chunk unnamed-chunk-4

Solution problème n°1

Affichage des résultats du modèle m1

summary(m1.1, prob = 0.89)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + gender 
   Data: d1 (Number of observations: 40) 
Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
         total post-warmup samples = 12000

Population-Level Effects: 
          Estimate Est.Error l-89% CI u-89% CI Eff.Sample Rhat
Intercept    63.80      0.70    62.70    64.90      11622 1.00
genderM       4.17      0.98     2.59     5.76      10153 1.00

Family Specific Parameters: 
      Estimate Est.Error l-89% CI u-89% CI Eff.Sample Rhat
sigma     3.11      0.37     2.57     3.74       9690 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

En Utilisant la library broom

library(broom)
tidy(m1.1, prob = 0.89)

         term    estimate std.error       lower       upper
1 b_Intercept   63.801407 0.6954508   62.696799   64.896623
2   b_genderM    4.172265 0.9794507    2.585206    5.761208
3       sigma    3.105483 0.3720946    2.574076    3.743628
4        lp__ -109.747336 1.2538574 -112.127615 -108.395307

En Utilisant la fonction posterior_summary

posterior_summary(m1.1)[1:3, ] %>%
    round(digits = 3)

            Estimate Est.Error   Q2.5  Q97.5
b_Intercept   63.801     0.695 62.429 65.174
b_genderM      4.172     0.979  2.249  6.116
sigma          3.105     0.372  2.488  3.943

Solution problème n°1

Comparaison des modèles

m1.1 <- add_criterion(m1.1, "waic")
m1.2 <- add_criterion(m1.2, "waic")
m1.3 <- add_criterion(m1.3, "waic")
m1.4 <- add_criterion(m1.4, "waic")

w <- loo_compare(m1.1, m1.2, m1.3, m1.4, criterion = "waic")

w[,5:8] %>% round(digits = 2)

     p_waic se_p_waic   waic se_waic
m1.3   5.17      1.47 187.01   10.67
m1.2   4.94      1.43 187.07   10.78
m1.4   5.33      1.57 188.06   10.90
m1.1   2.69      0.66 205.70    8.49

Poids respectifs des modèles

model_weights(m1.1, m1.2, m1.3, m1.4, weights = "waic") %>% 
  round(digits = 2)

m1.1 m1.2 m1.3 m1.4 
0.00 0.38 0.39 0.23

Solution problème n°1

my_coef_tab <-
  tibble(model = c("m1.1", "m1.2", "m1.3", "m1.4")) %>% 
  mutate(fit   = purrr::map(model, get)) %>% 
  mutate(tidy  = purrr::map(fit, tidy)) %>% 
  unnest(tidy) %>% 
  filter(term != "lp__")

my_coef_tab %>%
  complete(term = distinct(., term), model) %>%
  select(model, term, estimate) %>%
  mutate(estimate = round(estimate, digits = 2)) %>%
  spread(key = model, value = estimate)

# A tibble: 7 x 5
  term              m1.1   m1.2   m1.3   m1.4
  <chr>            <dbl>  <dbl>  <dbl>  <dbl>
1 b_father         NA     0.18    0.18  0.16 
2 b_genderM         4.17  3.54    7.3   1.3  
3 b_genderM:father NA    NA      NA     0.03 
4 b_genderM:mother NA    NA      -0.06 NA    
5 b_Intercept      63.8  15.1    13.7  16.1  
6 b_mother         NA     0.570   0.6   0.570
7 sigma             3.11  2.38    2.37  2.4

Solution problème n°1

# plot
my_coef_tab %>%
  ggplot(aes(x = model, y = estimate, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0,  color = "steelblue") +
  geom_pointrange(shape = 21, color = "steelblue", fill = "steelblue") +
  labs(x = NULL,
       y = NULL) +
  coord_flip() +
  theme_classic(base_size  = 30) +
  theme(text         = element_text(family = "Courier"),
        panel.grid   = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.y  = element_text(hjust = 0),
        panel.background = element_rect(),
        strip.background = element_rect(color = "transparent")) +
  facet_wrap(~term, ncol = 1)

plot of chunk unnamed-chunk-13

Solution problème n°1

pp_check(m1.3, nsamples = 1e2) + theme_bw(base_size = 20)

plot of chunk unnamed-chunk-14

Problème n°2

Les données suivantes documentent le naufrage du titanic. La colonne pclass indique la classe dans laquelle chaque passager voyageait (un proxy pour le statut socio-économique), tandis que la colonne parch indique le nombre de parents et enfants à bord.

Peut-on prédire la survie d'un passager en fonction de ces informations ?

d2 <- read.csv("Data/titanic.csv")
glimpse(d2)

Observations: 539
Variables: 5
$ survival <int> 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1…
$ pclass   <fct> upper, lower, upper, lower, upper, lower, upper, upper,…
$ gender   <fct> male, female, female, female, male, male, male, female,…
$ age      <dbl> 22, 38, 26, 35, 35, 54, 2, 27, 4, 58, 20, 39, 14, 2, 31…
$ parch    <int> 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 5, 0, 1, 0, 0, 0, 1, 5…

Solution problème n°2

Cette situation revient à essayer de prédire un outcome dichotomique à l'aide de prédicteurs continus et / ou catégoriels. On peut utiliser un modèle de régression logistique (cf. cours n°6).

# centering and scaling
d2 <-
    d2 %>%
    mutate(
        pclass = ifelse(pclass == "lower", -0.5, 0.5),
        gender = ifelse(gender == "female", -0.5, 0.5),
        age = scale(age) %>% as.numeric,
        parch = scale(parch) %>% as.numeric
        )

Solution problème n°2

d2 %>%
    group_by(pclass, gender) %>%
    summarise(p = mean(survival) ) %>%
    ggplot(aes(x = as.factor(pclass), y = p, fill = as.factor(gender) ) ) +
    geom_bar(position = position_dodge(0.5), stat = "identity") +
    theme_bw(base_size = 25) +
    xlab("class") + ylab("p(survival)")

plot of chunk unnamed-chunk-17

Solution problème n°2

m0 <- brm(
    survival ~ 1,
    family = binomial(link = "logit"),
    prior = prior(normal(0, 10), class = Intercept),
    data = d2,
    cores = parallel::detectCores()
    )

m1 <- brm(
    # using the dot is equivalent to say "all predictors" (all columns)
    survival ~ .,
    family = binomial(link = "logit"),
    prior = c(
        prior(normal(0, 10), class = Intercept),
        prior(normal(0, 10), class = b)
    ),
    data = d2,
    cores = parallel::detectCores()
    )

m2 <- update(m1,
    formula = survival ~ 1 + pclass + gender + pclass:gender)

m3 <- update(m1,
    formula = survival ~ 1 + pclass + gender + pclass:gender + age)

Solution problème n°2

m3 %>%
    plot(
        pars = "^b_",
        combo = c("hist", "trace"), widths = c(1, 2.5),
        theme = theme_bw(base_size = 20)
        )

plot of chunk unnamed-chunk-19

Solution problème n°2

summary(m3)

 Family: binomial 
  Links: mu = logit 
Formula: survival ~ pclass + gender + age + pclass:gender 
   Data: d2 (Number of observations: 539) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept         0.32      0.18    -0.01     0.71       1551 1.00
pclass           -2.97      0.40    -3.81    -2.24       1463 1.00
gender           -2.61      0.36    -3.39    -1.98       1560 1.00
age              -0.47      0.13    -0.74    -0.23       2861 1.00
pclass:gender     2.26      0.73     0.97     3.79       1381 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

En Utilisant la library broom

library(broom)
tidy(m3, prob = 0.89)

             term     estimate std.error         lower        upper
1     b_Intercept    0.3243055 0.1846382    0.04603331    0.6296130
2        b_pclass   -2.9662498 0.3971888   -3.61622317   -2.3670721
3        b_gender   -2.6135923 0.3640584   -3.24780592   -2.0731990
4           b_age   -0.4748002 0.1294763   -0.68627376   -0.2724301
5 b_pclass:gender    2.2636830 0.7257193    1.19155467    3.5041731
6            lp__ -269.8048698 1.5945055 -272.82650516 -267.8832857

En Utilisant la fonction posterior_summary

posterior_summary(m3)[1:3, ] %>%
    round(digits = 3)

            Estimate Est.Error   Q2.5  Q97.5
b_Intercept    0.324     0.185 -0.009  0.711
b_pclass      -2.966     0.397 -3.806 -2.244
b_gender      -2.614     0.364 -3.387 -1.978

Solution problème n°2

m0 <- add_criterion(m0, "waic")
m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")

w <- loo_compare(m0, m1, m2, m3, criterion = "waic")

w[,5:8] %>% round(digits = 2)

   p_waic se_p_waic   waic se_waic
m3   5.33      0.80 512.69   26.52
m1   5.05      0.45 520.99   25.91
m2   4.24      0.66 524.75   25.42
m0   0.97      0.02 717.98   10.99

Poids respectifs des modèles

model_weights(m0, m1, m2, m3, weights = "waic") %>% 
  round(digits = 2)

  m0   m1   m2   m3 
0.00 0.02 0.00 0.98

Solution problème n°2

pp_check(m3, nsamples = 1e2) + theme_bw(base_size = 25)

plot of chunk unnamed-chunk-25

Problème n°3

Ce jeu de données recense des informations sur le diamètre (colonne diam) de 80 pommes (chaque pomme étant identifiée par la colonne id), poussant sur 10 arbres différents (colonne tree). On a mesuré ce diamètre pendant 6 semaines successives (colonne time).

Que peut-on dire de la pousse de ces pommes, tout en considérant les structures de dépendance existant dans les données (chaque pomme poussait sur un arbre différent) ?

d3 <- read.csv("Data/apples.csv")
glimpse(d3)

Observations: 480
Variables: 5
$ tree  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ apple <int> 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 10, …
$ id    <int> 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 10, …
$ time  <int> 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2…
$ diam  <dbl> 2.90, 2.90, 2.90, 2.93, 2.94, 2.94, 2.86, 2.90, 2.93, 2.96…

Solution problème n°3

Cette situation revient à essayer de prédire un outcome dichotomique à l'aide de prédicteurs continus et / ou catégoriels. On peut utiliser un modèle multi-niveaux (ou modèle mixte, cf. cours n°8).

d3 <- d3 %>% filter(diam != 0) # removing missing data

d3 %>% # plotting
    ggplot(aes(x = time, y = diam, colour = as.factor(apple) ) ) +
    geom_point(show.legend = FALSE) +
    geom_line(show.legend = FALSE) +
    facet_wrap(~tree, ncol = 5) +
    theme_bw(base_size = 25)

plot of chunk unnamed-chunk-27

Solution problème n°3

On peut fitter plusieurs modèles avec brms et les comparer ensuite en utilisant le WAIC.

p1 <- c(
    prior(normal(0, 10), class = Intercept),
    prior(cauchy(0, 10), class = sigma)
    )

m1 <- brm(
    diam ~ 1,
    prior = p1,
    data = d3,
    cores = parallel::detectCores()
    )

p2 <- c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma)
    )

m2 <- brm(
    diam ~ 1 + time,
    prior = p2,
    data = d3,
    cores = parallel::detectCores()
    )

Solution problème n°3

p3 <- c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sd),
    prior(cauchy(0, 10), class = sigma)
    )

m3 <- brm(
    diam ~ 1 + time + (1 | tree),
    prior = p3,
    data = d3,
    cores = parallel::detectCores()
    )

p4 <- c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sd),
    prior(cauchy(0, 10), class = sigma),
    prior(lkj(2), class = cor)
    )

m4 <- brm(
    diam ~ 1 + time + (1 + time | tree),
    prior = p4,
    data = d3,
    cores = parallel::detectCores(),
    control = list(adapt_delta = 0.99)
    )

Solution problème n°3

p5 <- c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sd),
    prior(cauchy(0, 10), class = sigma),
    prior(lkj(2), class = cor)
    )

m5 <- brm(
    diam ~ 1 + time + (1 + time | tree / apple),
    prior = p5,
    data = d3,
    cores = parallel::detectCores(),
    control = list(adapt_delta = 0.99)
    )

Solution problème n°3

WAICs

m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")
m4 <- add_criterion(m4, "waic")
m5 <- add_criterion(m5, "waic")

w <- loo_compare(m1, m2, m3, m4, m5, criterion = "waic")

w[,5:8] %>% round(digits = 2)

   p_waic se_p_waic     waic se_waic
m5 114.76      7.92 -2297.34   33.76
m3  12.25      1.54  -824.42   43.24
m4  13.71      1.66  -822.92   43.26
m2   4.38      1.02  -778.63   48.58
m1   2.94      0.77  -701.90   43.54

Poids respectifs des modèles

model_weights(m1, m2, m3, m4, m5, weights = "waic") %>% 
  round(digits = 2)

m1 m2 m3 m4 m5 
 0  0  0  0  1

Solution problème n°3

summary(m5)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: diam ~ 1 + time + (1 + time | tree/apple) 
   Data: d3 (Number of observations: 451) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~tree (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)           0.02      0.01     0.00     0.05        339 1.00
sd(time)                0.00      0.00     0.00     0.01       1703 1.00
cor(Intercept,time)     0.20      0.42    -0.68     0.88        708 1.00

~tree:apple (Number of levels: 80) 
                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)           0.09      0.01     0.08     0.11        846 1.00
sd(time)                0.01      0.00     0.00     0.01       1630 1.00
cor(Intercept,time)     0.50      0.13     0.23     0.73       2201 1.00

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     2.83      0.01     2.81     2.86        661 1.01
time          0.03      0.00     0.02     0.03       1902 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.02      0.00     0.01     0.02       2624 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Solution problème n°3

post <- posterior_samples(m5, "b")

ggplot(data = d3, aes(x = time, y = diam) ) +
    geom_point(alpha = 0.5, shape = 1) +
    geom_abline(
        data = post, aes(intercept = b_Intercept, slope = b_time),
        alpha = 0.01, size = 0.5) +
    theme_bw(base_size = 25)

plot of chunk unnamed-chunk-34

Solution problème n°3

library(tidybayes)
library(modelr)

d3 %>%
    group_by(tree, apple) %>%
    data_grid(time = seq_range(time, n = 1e2) ) %>%
    add_fitted_samples(m5, n = 1e2) %>%
    ggplot(aes(x = time, y = diam, colour = factor(apple) ) ) +
    geom_line(
        aes(y = estimate, group = paste(apple, .iteration) ),
        alpha = 0.2, show.legend = FALSE) +
    facet_wrap(~tree, ncol = 5) +
    theme_bw(base_size = 25)

$plot of chunk unnamed-chunk-35$

Solution problème n°3

pp_check(m3, nsamples = 1e2) + theme_bw(base_size = 25)

plot of chunk unnamed-chunk-36

pp_check(m5, nsamples = 1e2) + theme_bw(base_size = 25)

plot of chunk unnamed-chunk-37

Problème n°4

Ces données recensent le nombre de candidatures pour 6 départements (colonne dept) à Berkeley (données disponibles dans le package rethinking). La colonne admit indique le nombre de candidatures acceptées et la colonne reject le nombre de candidatures rejetées (la colonne applications est simplement la somme des deux), en fonction du genre des candidats (applicant.gender).

On veut savoir s'il existe un biais de genre dans l'admission des étudiants à Berkeley.

d4 <- read.csv("Data/UCBadmit.csv")
glimpse(d4)

Observations: 12
Variables: 6
$ X                <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
$ dept             <fct> A, A, B, B, C, C, D, D, E, E, F, F
$ applicant.gender <fct> male, female, male, female, male, female, male,…
$ admit            <int> 512, 89, 353, 17, 120, 202, 138, 131, 53, 94, 2…
$ reject           <int> 313, 19, 207, 8, 205, 391, 279, 244, 138, 299, …
$ applications     <int> 825, 108, 560, 25, 325, 593, 417, 375, 191, 393…

Solution problème n°4

Cette situation revient à essayer de prédire un outcome dichotomique (admit, reject) à l'aide de prédicteurs continus et / ou catégoriels.

d4 %>%
    ggplot(aes(x = dept, y = admit / applications) ) +
    geom_bar(stat = "identity") +
    facet_wrap(~applicant.gender) +
    theme_bw(base_size = 25)

plot of chunk unnamed-chunk-39

Solution problème n°4

On peut fitter plusieurs modèles avec brms et les comparer ensuite en utilisant le WAIC.

# centering gender predictor
d4$gender <- ifelse(d4$applicant.gender == "female", -0.5, 0.5)

# creating an index for department
d4$dept_id <- coerce_index(d4$dept)

p1 <- c(
    prior(normal(0, 10), class = "Intercept"),
    prior(cauchy(0, 2), class = "sd")
    )

m1 <- brm(
    admit | trials(applications) ~ 1 + (1 | dept_id),
    data = d4, family = binomial,
    prior = p1,
    warmup = 1000, iter = 5000,
    control = list(adapt_delta = 0.99, max_treedepth = 12)
    )

p2 <- c(
    prior(normal(0, 10), class = "Intercept"),
    prior(normal(0, 1), class = "b"),
    prior(cauchy(0, 2), class = "sd")
    )

m2 <- brm(
    admit | trials(applications) ~ 1 + gender + (1 | dept_id),
    data = d4, family = binomial,
    prior = p2,
    warmup = 1000, iter = 5000,
    control = list(adapt_delta = 0.99, max_treedepth = 12)
    )

Solution problème n°4

p3 <- c(
    prior(normal(0, 10), class = "Intercept"),
    prior(normal(0, 1), class = "b"),
    prior(cauchy(0, 2), class = "sd"),
    prior(lkj(2), class = "cor")
    )

m3 <- brm(
    admit | trials(applications) ~ 1 + gender + (1 + gender | dept_id),
    data = d4, family = binomial,
    prior = p3,
    warmup = 1000, iter = 5000,
    control = list(adapt_delta = 0.99, max_treedepth = 12)
    )

Solution problème n°4

m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")

w <- loo_compare(m1, m2, m3, criterion = "waic")

w[,5:8] %>% round(digits = 2)

   p_waic se_p_waic   waic se_waic
m3   6.85      1.32  91.15    4.60
m1   6.51      2.28 105.13   17.95
m2   9.37      3.04 108.39   16.54

model_weights(m1, m2, m3, weights = "waic") %>%
  round(digits = 2)

m1 m2 m3 
 0  0  1

stanplot(m3, pars = c("^b_", "^sd_")) +
  theme_fivethirtyeight(base_size = 25) +
  theme(axis.text.y = element_text(hjust = 0))

plot of chunk unnamed-chunk-43

library(broom)
tidy(m3, par_type = "non-varying", prob = 0.95)

       term   estimate std.error      lower     upper
1 Intercept -0.5800742 0.6772113 -1.9556560 0.7947479
2    gender -0.1622283 0.2471726 -0.6639525 0.3396001

Solution problème n°4

library(tidybayes)
library(modelr)

d4 %>%
    group_by(dept_id, applications) %>%
    data_grid(gender = seq_range(gender, n = 1e2) ) %>%
    add_fitted_samples(m3, newdata = ., n = 100, scale = "linear") %>%
    mutate(estimate = plogis(estimate) ) %>%
    ggplot(aes(x = gender, y = estimate, group = .iteration) ) +
    geom_hline(yintercept = 0.5, lty = 2) +
    geom_line(aes(y = estimate, group = .iteration), size = 0.5, alpha = 0.1) +
    facet_wrap(~dept_id, nrow = 2) +
    theme_bw(base_size = 20)

plot of chunk unnamed-chunk-45

Solution problème n°4

pp_check(m3, nsamples = 1e2) + theme_bw(base_size = 25)

plot of chunk unnamed-chunk-46

Problème n°5

Le dilemme du tramway (trolley problem) est une expérience de pensée qui permet d'étudier les déterminants des jugements de moralité (i.e., qu'est-ce qui fait qu'on juge une action comme morale, ou pas ?).

Sous une forme générale, ce dilemme consiste à poser la question suivante: si une personne peut effectuer un geste qui bénéficiera à un groupe de personnes A, mais, ce faisant, nuira à une personne B (seule); est-il moral pour la personne d'effectuer ce geste ?

Voir ce lien pour plus d'informations.

plot of chunk unnamed-chunk-47

Problème n°5

Généralement, on fait lire des scénarios aux participants de l'étude, dans lesquels un individu doit prendre une décision dans une situation similaire à celle décrite à la slide précédente. Par exemple, imaginons que Denis ait le choix entre ne rien faire et laisser un train tuer cinq personnes, ou faire dérailler ce train mais tuer une personne. Ensuite, on demande aux participants de juger de la moralité de l'action choisie de Denis, sur une échelle de 1 à 7.

Des études antérieures ont montré que ces jugements de moralité sont grandement influencés par trois mécanismes de raisonnement inconscients:

Le principe d'action: un préjudice causé par une action est jugé moralement moins acceptable qu'un préjudice causé par omission
Le principe d'intention: un préjudice causé come étant le moyen vers un but est jugé moralement moins acceptable qu'un préjudice étant un effet secondaire (non désiré) d'un but
Le principe de contact: un préjudice causé via contact physique est jugé moralement moins acceptable qu'un préjudice causé sans contact physique

Problème n°5

Ce jeu de données comprend 12 colonnes et 9930 lignes, pour 331 individus. L'outcome qui nous intéresse est response, un entier pouvant aller de 1 à 7, qui indique à quel point il est permis (moralement) de réaliser l'action décrite dans le scénario correspondant, en fonction de l'age et genre (male) du participant (id).

On se demande comment les jugements d'acceptabilité sont influencés par les trois principes décrits slide précédente. Ces trois principes correspondent aux trois variables, action, intention, et contact (dummy-coded).

d5 <- read.csv("Data/morale.csv")
glimpse(d5)

Observations: 9,930
Variables: 7
$ response  <int> 4, 3, 4, 3, 3, 3, 5, 4, 4, 4, 4, 4, 4, 5, 4, 4, 4, 4, …
$ id        <fct> 96;434, 96;434, 96;434, 96;434, 96;434, 96;434, 96;434…
$ age       <int> 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14…
$ male      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ action    <int> 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, …
$ intention <int> 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ contact   <int> 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …

Solution problème n°5

On essaye de prédire un jugement exprimé sous forme d'entier entre 1 et 7. Autrement dit, la variable qu'on essaye de prédire est une variable catégorielle, et dont les catégories sont ordonnées de 1 à 7…

d5$response %>% table %>% plot(xlab = "response", ylab = "", cex.axis = 2, cex.lab = 2)

plot of chunk unnamed-chunk-49

Solution problème n°5

Ce type de données peut se modéliser en utilisant le modèle de régression logistique ordinale. Ci-dessous un exemple en utilisant brms, et les priors par défaut (NB: ces modèles peuvent être un peu longs à fitter…).

moral1 <- brm(
    response ~ 1,
    data = d5,
    family = cumulative("logit"),
    cores = parallel::detectCores(),
    )

moral2 <- brm(
    response ~ 1 + action + intention + contact,
    data = d5,
    family = cumulative("logit"),
    cores = parallel::detectCores()
    )

Solution problème n°5

Toutes les pentes sont négatives… ce qui signifie que chaque facteur réduit la réponse moyenne (i.e., le jugement de moralité). Ces pentes représentent des changements dans les log-odds cumulatifs.

moral1 <- add_criterion(moral1, "waic")
moral2 <- add_criterion(moral2, "waic")

w <- loo_compare(moral1, moral2, criterion = "waic")

w[,5:8] %>% round(digits = 2)

       p_waic se_p_waic     waic se_waic
moral2   9.13      0.04 37090.12   76.28
moral1   6.01      0.02 37854.46   57.67

model_weights(moral1, moral2, weights = "waic") %>%
  round(digits = 2)

moral1 moral2 
     0      1

tidy(moral2, prob = 0.95)

             term      estimate  std.error         lower         upper
1  b_Intercept[1] -2.839163e+00 0.04805513 -2.932892e+00 -2.744327e+00
2  b_Intercept[2] -2.156854e+00 0.04278314 -2.240424e+00 -2.074097e+00
3  b_Intercept[3] -1.574014e+00 0.03950732 -1.650146e+00 -1.498460e+00
4  b_Intercept[4] -5.524352e-01 0.03699212 -6.240745e-01 -4.799432e-01
5  b_Intercept[5]  1.164427e-01 0.03713198  4.442318e-02  1.895871e-01
6  b_Intercept[6]  1.023401e+00 0.03987168  9.436735e-01  1.100674e+00
7        b_action -7.093860e-01 0.04046395 -7.861238e-01 -6.308170e-01
8     b_intention -7.212895e-01 0.03688297 -7.940994e-01 -6.505540e-01
9       b_contact -9.614392e-01 0.05087764 -1.063520e+00 -8.625872e-01
10           lp__ -1.856178e+04 2.11921873 -1.856673e+04 -1.855863e+04

Solution problème n°5

On peut représenter les prédictions du modèle en utilisant la fonction marginal_effects.

marg1 <- marginal_effects(moral2, "action", ordinal = TRUE)

p1 <-
    plot(marg1, theme = theme_bw(base_size = 25), plot = FALSE)[[1]] +
    scale_y_discrete(breaks = 1:7, labels = 1:7)

marg2 <- marginal_effects(moral2, "intention", ordinal = TRUE)

p2 <-
    plot(marg2, theme = theme_bw(base_size = 25), plot = FALSE)[[1]] +
    scale_y_discrete(breaks = 1:7, labels = 1:7)

marg3 <- marginal_effects(moral2, "contact", ordinal = TRUE)

p3 <-
    plot(marg3, theme = theme_bw(base_size = 25), plot = FALSE)[[1]] +
    scale_y_discrete(breaks = 1:7, labels = 1:7)

Solution problème n°5

ggpubr::ggarrange(p1, p2, p3, nrow = 1, ncol = 3, common.legend = TRUE)

plot of chunk unnamed-chunk-53

Solution problème n°5

Posterior predictive checking. Pour plus d'informations sur la régression logistique ordinale, voir ce lien, ce tutorial, ou le chapitre 11 de McElreath (2015).

pp_check(moral2, nsamples = 1e2) + theme_bw(base_size = 25)

plot of chunk unnamed-chunk-54